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Abstract. We discuss the preparation of many-body states of cold fermionic atoms in 
an optical lattice via controlled dissipative processes induced by coupling the system to 
a reservoir. Based on a mechanism combining Pauli blocking and phase locking between 
adjacent sites, we construct complete sets of jump operators describing coupling to 
a reservoir that leads to dissipative preparation of pairing states for fermions with 
various symmetries in the absence of direct inter-particle interactions. We discuss the 
uniqueness of these states, and demonstrate it with small-scale numerical simulations. 
In the late time dissipative dynamics, we identify a "dissipative gap" that persists 
in the thermodynamic limit. This gap implies exponential convergence of all many- 
body observables to their steady state values. We then investigate how these pairing 
states can be used as a starting point for the preparation of the ground state of Fermi- 
Hubbard Hamiltonian via an adiabatic state preparation process also involving the 
parent Hamiltonian of the pairing state. We also provide a proof-of-principle example 
for implementing these dissipative processes and the parent Hamiltonians of the pairing 
states, based on ^''^Yb atoms in optical lattice potentials. 
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1. Introduction 

Quantum simulation using cold atoms in an optical lattice typically requires cooling 
to low temperatures to see interesting quantum phases with strong correlations 
[H El El SI IS]- An important example in this context is the quantum simulation of 
the 2-D Fermi-Hubbard Model (FHM) using cold fermionic atoms in an optical lattice 
[6]. As a crucial first step toward this goal, various experimental groups have been 
successful in realising the 3-D Fermi-Hubbard model in such a system [71 El El ITO] . 
But cooling of the system below the critical temperature, and thus into the phases of 
interest, turns out to be difficult with the conventional cooling schemes (T/Tp ~ 0.25 
[HI [To], where Tp is the Fermi temperature). On the other hand, the atomic physics 
of these systems opens the possibility for a different approach to the production of 
interesting many-body states, specifically to dissipatively drive the system into steady 
states with the desired coherence and symmetry properties by careful engineering of a 
reservoir [TH [T21 [131 HH ESI HSl 113 HSl US] ■ In this paper, we will discuss schemes by 
which fermions in a 2-D lattice potential can be dissipatively driven into pairing states 
with non-trivial correlations even in the absence of attractive interactions, extending 
our previous work [T9]. These many-body pairing states can then also be used as 
a low-entropy starting point for efficient adiabatic passages, through which the true 
many-body ground state of many-body models such as the Fermi-Hubbard model might 
be reached [31 [115]. 

For Markovian dissipative processes, the coupling with the reservoir can be modeled 
by a set of dissipative quantum jump operators, which, when chosen appropriately, drive 
the system towards a pure steady state with zero entropy starting from any given initial 
state. Similar ideas have been applied recently for the dissipative preparation of many- 
body states in optical lattices, e.g. a Bose-Einstein condensate for bosons [12], an rj- 
condensate for fermions [13] , and most recently, d-wave pairing states in optical lattices 
[19] as well as topological phases [20] . Here, we will focus on the dissipative generation of 
pairing states for fermions in an optical lattice. Starting from the general principles, we 
will show in detail how relative phases between the atoms can be imposed by engineering 
the jump operators, which in turn allows us to prepare many-body pairing states with 
specific spatial symmetries, e.g. p-wave and d-wave pairing symmetries. Note that in 
contrast to the equilibrium states of typical non-dissipative processes, the dissipative 
dynamics is described by a master equation, and the final state that we aim to prepare 
appears as the steady state of the dynamical process, providing a targeted many-body 
cooling protocol. As a conceptually important result, we observe that the imaginary 
spectrum of the effective Hamiltonian of the master equation has a gap which persists 
in the thermodynamic limit and we refer to it as a "dissipative gap" , due to the formal 
analogy to the pairing gap of conventional BCS paired states [2T]. Physically, it leads 
to exponential convergence of many-body observables to their steady state values in the 
late time dissipative dynamics. This is a unique feature of the dissipative preparation 
of paired states with fermions. 
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We then propose an adiabatic process to connect the steady state of the dissipative 
process with the ground state of a Hamiltonian with similar symmetry properties. In 
the case of ideal adiabaticity, the system remains in a pure state and evolves into the 
ground state of the Hamiltonian without going through any Landau-Zener crossings. 
We consider this in the context of the Fermi-Hubbard model, where a dissipative 
process can be used to create dissipatively bound fermion pairs which have d-wave 
symmetry. As we will see later, this BCS-type mean-field d-wave state lacks the 
strong correlations associated with the repulsive Fermi-Hubbard model, as the double- 
occupancy is not projected out by the dissipative dynamics. We will show that strong 
correlations can be built up through proper adiabatic passage in the sense that the 
steady state can be connected with the ground state of the Fermi-Hubbard model 
efficiently. The central question in the study of the 2-D FHM is whether the ground 
state exhibits d-wave superconductivity /superfluidity away from half-flUing, and thus 
captures the universal properties shared by the high-Tc superconducting materials 
[SaESlElESlEniETlEHlEnlEO]. if this is indeed the case, the coherence and the 
d-wave symmetry should be conserved during the adiabatic process. 

Finally, we note that a key physical ingredient for the jump operators for pairing 
states of fermions is the Pauli exclusion principle [31]. Based on this understanding, 
we propose to implement the jump operators stroboscopically using alkaline-earth-like 
atoms. The existence of long-lived meta-stable states and rich level structures in these 
atoms provides us with the freedom to engineer various dissipative processes. 

The paper is organised as follows: In Sec. II, we discuss the general formalism 
for jump operator engineering, and derive both the fixed-number and the fixed-phase 
jump operators for pairing states with different symmetries; in Sec. HI, we study the 
uniqueness of the steady state under these jump operators for various cases both from 
the symmetry perspective and with small scale numerical calculations; we then derive a 
mean-field theory in Sec. IV for the master equation describing the dissipative dynamics, 
where a "dissipative gap" for the pairing states - a minimal damping rate for the many- 
body relaxation - is shown to emerge; in Sec. V, having established the dissipative 
preparation of pure steady states, we illustrate how we may start from these initial 
states to prepare the ground state of Hamiltonians with similar symmetries, for example 
the FHM, via an efficient adiabatic process; we discuss the implementation of the jump 
operators using alkaline-earth-like atoms in Sec. VI; finally we summarise and discuss 
other possible apphcations of our dissipative state preparation setting. 

2. General principles of reservoir engineering 

Closed systems, where the energy and particle number are conserved, are modeled by a 
Hamiltonian, which determines the ground state and the dynamics via the Schrodinger 
equation. In this context, engineering states by realisation of a particular Hamiltonian, 
for which the desired state is the ground state is often discussed. It is then natural to 
ask similar questions for open quantum systems, where dissipative processes interrupt 
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coherent evolution. In particular, one can consider engineering the dissipative processes 
so that their back-actions project the system into the desired subspace of the complete 
Hilbert space. For Markovian dissipative processes, as appropriate for the systems 
considered below, the dynamics of the density matrix for an open system is described 
by a master equation. 



Here, {je} are non-Hermitian Lindblad operators reflecting the system-bath coupling 
with rate k [3^. The Hamiltonian H generates unitary evolution, and describes non- 
dissipative processes of the system. Although H does not have to vanish in general, we 
will assume H = throughout the state preparation discussion, i.e. the final state is 
prepared via purely dissipative processes. For fermions loaded into an optical lattice, 
this can be achieved by increasing the lattice depth to freeze out the kinetic motion, 
while tuning the inter-particle interaction to zero, e.g. via a Feshbach resonance. Note 
that by considering a purely dissipative process, we can avoid competition between the 
Hamiltonian and the dissipative dynamics, which is present when the dark state is not 
an exact eignestate of the Hamiltonian. In particular, the pairing states that we aim to 
prepare are not exact eigenstates of the FHM in general. 

In the quantum trajectory picture, the system wavefunction {ipit)) of a given 
trajectory evolves according to the non-Hermitian Hamiltonian lipit)) oc e~''^''^^\xl>{0)), 
and is punctuated by the quantum jump \ip{t)) — )■ je\'ip{t)) with rate k, \\je\4'{t))\f ■ The 
time-dependent density matrix is then determined by p(t) = {\4'(t)){ip(t)\)stoch [S3], 
where the average runs over all trajectories. In this picture, we see that for any state 
satisfying j^lBCSiv) = V£, the quantum jumps will never project it to other states. 
These states are therefore "dark states" of the jump operators, and are necessarily 
steady states of the master equation evolution [121 E] • If other stationary solutions 
exist, the system will be driven to this state by the dissipative dynamics regardless of its 
initial conditions [11]. Similar ideas have been exploited for the dynamical preparation 
of a BEC for cold bosons loaded into an optical lattice, where the jump operators give 
rise to quasi-local phase-locking mechanism, which eventually leads to the condensation 
of the bosons in the lattice [10]. Here, we will focus on the preparation of pairing 
states for fermions in an optical lattice. For the fermionic case here, in addition to the 
phase-locking mechanism as in the case of bosons, the physical foundation of the jump 
operators is the Pauli blocking (see Fig. [T]), i.e. for fermions the spontaneous emission 
to an already occupied state is blocked [31]. This gives rise to a novel non-equilibrium 
pairing mechanism for fermions that does not require attractive conservative forces. 




(1) 



where the non-Hermitian effective Hamiltonian is given by 




(2) 
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Figure 1. Illustration of the Pauli blocking in the optical pumping process of a three 
level system, (a) The spontaneous decay from the excited state |e) to the ground state 
1^2) is blocked due to Pauli exclusion, leaving the system unchanged; (b) The decay 
channel is not blocked, and the population in \gi) is transferred to the state \g2)- 

We are primarily interested in states with a homogeneous product of identical 
fermion pairs on a two-dimensional (2-D) square lattice: 

|^)=r/t^|vac), = J2va,bclcl (3) 

a,b 

Here, c], creates a fermion in mode a, where a = (a, i) labels spin a and position i on 
the 2-D lattice. Typically, when we consider delocalised states, the sum extends over 
the whole lattice in the position space index. For large systems in the thermodynamic 
limit, we may adopt the grand canonical ensemble, and the state above becomes the 
BCS-type coherent state for paired fermions: 

|v]>) =Arexp(^/^,Bcl44)|vac) =ArJ](l + /^,Bcli4)|vac), (4) 

A,B A,B 

where we have Fourier transformed Eq. (|3| into momentum space, with A = {a, k} 
labeling spin a and momentum k, Af being the normalization factor. The pair 
wavefunction in momentum space is fA,B = jVa,b^^^''''^^^''^^ , where the spin degrees 
of freedom are not changed. The symmetry of the pairing state is encoded in the form 
of the pair wavefunction rja^ {fA,B)- 

In general, the jump operator that drives the system into the pairing state of Eq. 
Q can involve multi-particle processes. Single-particle jump operators (involving at 
most one annihilation mode operator), which act directly on one particle at a time, 
are generally easier to implement experimentally than two-particle jump operators 
(involving two annihilation mode operators). However, the two-particle jump operators 
are the most intuitive, as the dissipative dynamics can be viewed as local phase- 
locking between pairs of fermions. We will therefore first give a brief description on 
the derivation of two-particle jump operators for the pairing states, which will provide 
important clues as to how to proceed with the more practical single-particle operators. 
In both cases, our focus will be on the d-wave pairing state, though the procedures can 
be extended to other symmetries as well. 

The fixed-number d-wave pairing state is given by the symmetric superposition of 
d-wave pairs, which are spin-singlet fermion pairs on the bonds in x and y direction 
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Figure 2. Schematics for tlie symmetry of antiferromagnetic Neel state and d-wave 
state, (a) Spin configurations for Neel state. Particles with different colors represent 
different spins; (b) d-wave symmetry on a lattice. Each leaf for the clover-leaf structure 
here represents singlet pairing of spins on the adjacent sites. The singlet pairs change 
sign when rotated by 7r/4, as dictated by the d-wave symmetry. 



whose relative phase changes sign under a 7r/4 rotation x ^ y, 
= B»"|vac), 



(5) 
(6) 



with e^;, Cy the unit lattice vectors long the x and y direction, respectively. 

Following Eq. Q , it is easy to write down the general form of a fixed-phase d-wave 
pairing state: 



A/ > — I vac 



n 



) = J\f exp (aD^) |vac), 



(7) 



where a = e*^|a| is a complex number carrying the phase 6 of the pairing state. In 
the following, we will first focus on the fixed-number state and will come back to the 
fixed-phase state later. 

Assuming translational invariance, which is the case for an infinitely large 
homogeneous lattice or for a finite homogeneous lattice with periodic boundary 
conditions, we may rewrite the pair creation operators: 

4 = {(sVe.t + S-e.t) - (sVe,t + (8) 

SO that Z^^ = dj. These operators have the advantage that they are already factorised 
in the spin degrees of freedom, which allows us to write the pairing jump operators in 
a similar fashion as well. For convenience, we adopt a shorthand convention and write 
= J2u Pi'(^]+e^^(^j,iy where p±x = 1, p±y = —1. We also find it useful to define the 
singlet pairing state in one dimension (1-D): 



(r/t)^|vac) 



-I N 



El 



I vac). 



(9) 



where c||^^| is the single fermion annihilation operator on the ith site with a certain 
spin. The 1-D state [i]^)'^ |vac) captures the off-site singlet pairing feature of the d-wave 
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state. Similar to the d-wave case in 2-D, we may use translational symmetry to define 
the 1-D singlet pairing operator rjj: 

so that f]^ = J2jVj- 

In the following, we proceed by devising the jump operators for the 1-D singlet 
pairing state first before considering the 2-D scenarios where additional spatial phase 
locking is needed. 



2.1. Two-particle jump operators 

We are interested in jump operators that conserve particle number. This implies that 
they must carry total charge (or with total global phase 0). Thus, we can write it in 
normal order as a product of a pure creation and a pure annihilation part. The most 
general form of the two-particle jump operator thus reads 

Ji = xki, Xl = ^X*a,b(^)cicl ^i = ^Ub(j')CaCb. (H) 

a,b a,b 

Here we also impose the requirement of quasi-locality, i.e. the functions Xabi'i) , ^ab{i) 
shall be non-zero only in a small vicinity of site i in position space. Thus the jump 
operator Jj is centered around site i. 

To uniquely drive the system into the desired pairing state, it is necessary that the 
state |\E') be a dark state of a set of jump operators {xl^i} 

J,(r/t)^|vac) = xlUvYl^^c) = 0, Vz. (12) 

For a given operator ^j, we can work out its commutation relation with the creation 
operator of the pairs: 

[e„r/t] = A, [A,v^]=B. (13) 

While A carries charge and is composed of a constant plus a normal-ordered 
second order term, i? is a superposition of pair creation operators and carries charge 2, 
which implies that [-B,?]^] = (cf. Appendix A). 

With these relations, we find that the commutator of with the homogeneous 
product (?7^) is characterised by the commutators A and B only: 

UvY = ivY^^ + Nir^r-'A + ^^^^^Bir^Y-'- (M) 

The first term on the right-hand side of the equation above gives zero when acting on the 
vacuum. Similarly, the normal-ordered part of A yields zero on the vacuum. Therefore, 
in order to satisfy the dark state condition, we need to find a set of quasi-local bilinear 
operators and xl which for a given rj'^ uniquely solve the two equations 

A|vac) = 0, xIb = 0. (15) 



For the 1-D case, after some derivation following the arguments leading to Eq. (|15|) 

(see 



Appendix A for details), we find a two-particle jump operator of the form: 

J, = C^MQt, M = (4,^ + c,t_iJ(Q+^ - (16) 
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where is an arbitrary superposition of single-fermion creation operators with spin-up. 
Note that as a spin flip operation only changes the overall sign, jump operators similar 
to Ji but with spins flipped also have the singlet pairing state Eq. ([o]) as a dark state. 
This is also true for the jump operators of the pairing states that we consider in this 
work, as similar symmetries in the spin degrees of freedom hold for all the pairing states 
that we will consider. 



For the 2-D case, we find the jump operator (see Appendix A for details) 



V 

where p±x = 1, P±y = —1- This operator has an interesting structure. It may be seen 
as a conditional dissipative process. The annihilation of a fermion sitting at site i must 
only take place if a superposition of fermions located at sites i + e^.i + Cy is coherently 
transferred to a superposition on the four sites {i ± e^,., i ± e^^} centered around site i. 

The two-particle jump operators above have the desired pairing states as the dark 
state. However, numerical simulation shows that in most cases they do not guarantee 
a unique dark state. Furthermore, these jump operators involve correlated dissipation 
of two particles and are therefore difficult to implement experimentally. Nevertheless 
the construction scheme above provides clues for the design of jump operators that are 
easier to implement and with a unique dark state. 

2. 2. Single-particle jump operators 

In this subsection, we will show that counter-intuitively, it is possible to design a set 
of jump operators for which only single particle operations are required. The single- 
particle jump operators dissipatively drive fermions into pairs as well as phase lock the 
pairs into the desired symmetry. More importantly, we will give arguments later that 
the dark state of these single-particle jump operators should be unique, and they are 
easier to implement than the two-particle jump operators. Below, we will describe two 
ways in which the appropriate single-particle jump operators can be derived. 

Firstly, consistent with the discussion in the previous section, in the case of a single- 
particle jump operator, the annihilation part contains a single annihilation operator 
which must be either Cj^i- or q^^.. Hence the operator A carries charge 1, and we have 
B = [A,r]^ = 0. Taking the 1-D singlet pairing state as an example, it is easy to derive 
that A = cl^-j^^^ + cj_i^4^ for = Ci^^; and A = cj^i^-^ + for = Thus to satisfy 

Eq. (12), we need = 0, and the simplest choice is xl = Ai+i,i + = '^iA'^ 



and Xi = (c|_,_i ^ + c|_^ ^) for = Ci^i- Considering the superpositions of these operators, 
we have the set of jump operators: 

Jf = (4i + 4_i)a'^c„ (18) 

with two-spinor Cj = (cj^i-, Cj,4.)^ and Pauli matrices cr" with a = ±,z. 

Alternatively, we start with the physical intuition that d-wave pairing states may be 
viewed as delocalised antiferromagnetic order away from half-filling. We may consider a 
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unit cell of the Neel state Sf^jlvac), where the Neel state unit cell operator Sf^j = cja'^Cj, 
and j is one of the nearest neighbors of site i which thus creates an adjacent pair 
of fermions with opposite spin. We notice that the singlet pairing operator in 1- 
D can be viewed as the superposition of two antiferromagnetic unit cell operators 
ril = S^j^_^_-^ + 5'/'j_]^. It is thus instructive to first construct the jump operators for 
an antiferromagnetic Neel state. To annihilate this unit cell Neel state, we simply need 
the Lindblad operators of the form: 

= cla^c,. (19) 

This jump operator generates hopping with a spin flip, which is impossible in the case 
that antiferromagnetic order is already present, due to the Fermi statistics. Generalising 
the unit cell operator to a 2-D lattice, we notice the Neel state can be written in 
eight different forms, |N±) = Hie^ "^J+eJ^ac) = (-l)^^/^ -q.^^ 5?^_^Jvac), with M the 
lattice size, and e^, = {ie^., ie^}. The 2-D Lindblad operators corresponding to those 
in Eq. (|l9|): 

jl,_,,^=cl,yc,,teAoTB. (20) 

These operators impose the quasi-local constraint on the steady state that any given 
site should have opposite spin with its nearest neighboring sites, which guarantees 
antiferromagnetic order at half-filling. However, there is still a two-fold degeneracy 
where the antiferromagnetic order differs by a total spin flip. However, there is still a 
two-fold degeneracy where the antiferromagnetic order differs by a total spin flip. As 
the Neel states cannot be reached by implementing the jump operators in Eq. () on any 
other states, and as the complete set of jump operators is invariant under a Hermitian 
conjugation, the dark subspace, (containing the Neel states) is isolated from the rest 
of the Hilbert space under the dissipative dynamics given by Eq. (). These problems 
can be solved, in principle, by adding a single jump operator, or more naturally a set of 
jump operators of the form: 

ji = 4,aCj,^, (21) 

for arbitrary i and its nearest neighbor j, and for either spin a =t,i- As an example, 
for 2 G A and a =1, the dark state is |N-|-). Note that comparing the jump operators 
in Eq. (19) with the unit cell operators 5*^^^-, we find that the jump operators can be 
obtained from the state generating unit cell operator via a particle-hole transformation 
c] — >■ Cj on the central site j. 

Now that we have found the Lindblad operators for the Neel state, we can proceed 
to generalise the operators to the d-wave case. First, we identify the "d-wave unit cell 
operators" : 

A— E/^^Ve.' (22) 

V 

where p±x = ~P±y = 1 and a = ±. We then perform a particle-hole transformation on 
the central site as in the case of the Neel state, and find operators: 

•Jj — Pi^Jti+e^, ' '^i ~ PuJli+e^ ■ (23) 
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It is easy to verify that [Jt^J2j^'j] = {a = ±, z), which is dictated by the Fermi 
statistics. Note that the d-wave coherence is estabhshed via quasi-local phase locking 
between adjacent cloverleaves of sites, as illustrated in Fig. |2} In this respect, the jump 
operators act similarly to those that establish phase coherence in bosonic systems [12]. 
For fermionic pairing states, however, as explained above, Pauli blocking is an additional 
key ingredient. 

We now consider the general case of quasi-local pairing states that can be factorised 
in the spin degrees of freedom. The unit cell operator of these pairing states can be 
expressed as superpositions of the Neel state unit cell operators: 

/3l = 5^P.cl+e.,.,c^, (24) 

V 

where the coefficients py encode the spatial symmetry of the pairs, and a\ and 
(T2 can be arbitrary combinations of spin configurations, including spinless fermions 
(Ti = 02- Performing the particle-hole transformation as above, we get the following 
jump operators: 

~ Pt^cl+e^,g-i'^»,g-2; (25) 

V 

which gives [Jj, = 0, so long as the coefficients satisfy the following equations: 

Xl^MP<^c]+e,,<xiCj+e,„ai = 0, for (Ti ^ ds, (26) 

XI P/.Pr.(cj+e. - c]-ejcj+e^ = 0' ^1 = ^2, (27) 

where the spin degrees of freedom have been factored out, and we have omitted the spin 



index in Eq. (27), as there is only one spin species. 



For the spinful case, it is easy to verify that Eq. (26) always holds regardless of 



the structure of p^. The spinless case though, is non-trivial in general. A particularly 



interesting example following Eq. (27) above is a 2-D p-wave state of spinless fermions 
generated by p^ ~ ^ivP^A-^-eyA with p^, = —p-x = —^Py = ip-y = 1, the Lindblad 
operators are Puc\^e^Ci} . The p-wave pairing state in 2-D can be prepared in two 
different chiralities (see Fig. ): Px + ipy and px — ipy, which shares the spatial symmetry 
with the pairing states in topological superconductors. However, as we will show later 
in the mean-field analysis, the p-wave states prepared in this way are still in the strong 
pairing limit, such that they are topologically trivial. The generation of topological order 
in 2-D systems will be discussed in a forthcoming publication. For stable dissipatively 
induced topological order in 1-D, see [20] . 



2.3. Jump operators for fixed-phase state 

In the previous discussion, we have focused on the jump operators for dark states in the 
form of Eq. ([s]), i.e. states with fixed total particle number. The dissipative processes 
characterised by these jump operators necessarily conserve the total particle number of 
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Px-iPy 



Figure 3. Schematics for the p-wave symmetry on a lattice for spinless fermions. Each 
leaf for the clover- leaf structure stands for triplet pairing of particles on the adjacent 
sites. The overall phase of the triplet pairs change by i or —i under a 7r/4 rotation, as 
dictated by the p-wave symmetry. 



the system, as is clear from the commutation relation [Jf A^] = 0, where is the total 
particle operator. The total particle number of the dark state in this case is given by 
that of the initial state. 



On the other hand, one can show (c.f. Appendix B) that for any given number- 



conserving pairing state which is a dark state of single-particle jump operators, one can 
always construct a set of linear jump operators that have the fixed-phase state in the 
form of Eq. ([7|) as a dark state. In the case of the d-wave pairing state, we consider the 
following jump operators 



3i,i 



PuC. 



Qc, 



(28) 
(29) 



where p±x = ~P±y = 1 as given by the d-wave symmetry, and P and Q are complex 
numbers. It is easy to show that jj,o-|^)d = 0, provided that P/Q = a, with the fixed- 
phase pairing state \il))d defined in Eq. ([T]). Note that the jump operators for the 
fixed-phase state do not conserve the total particle number, while the average particle 
density in the dark state is given as A^ = Xlq,^ (^ICq.aCq.alV')^ = Eqi^^M^' wheie 
V9(q) is the pair wavefunction in momentum space, and the summation over q runs over 
the first Brillouin zone. Importantly, the average particle density of the dark state here 
is determined by the parameters of the dissipative process, i.e. |a| can be chosen to fix 
a desired average density. 

It is straightforward to extend the analysis to pairing states with other symmetries. 
For the spinless p-wave pairing state for example, the jump operators for the fixed-phase 
dark state \ip)p = A/'exp(Q;p^)|vac) are 



3i 



P'^Pucl+^^ + Qci, 



(30) 



with P/Q = 2a, and the factor of 2 is due to the triplet pairing symmetry of the p-wave 
pairing state. For a more detailed discussion on jump operators for fixed-phase states. 
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see Appendix B 



3. Completeness of the jump operators 

To achieve the dissipative preparation of many-body states, the final steady state of 
the master equation should be unique. This requires that (i) the dark state be unique; 
(ii) the non-existence of stationary solutions other than this dark state [13]. Taking 
the 1-D singlet pairing state as an example, we shall first analyse the uniqueness of the 
dark state from the symmetry perspective. We will then illustrate the uniqueness of the 
steady state with different pairing symmetries both in 1-D and 2-D by directly evolving 
the master equation for small finite size systems. 

3.1. Uniqueness of dark state from a symmetry perspective 

The problem of showing the uniqueness for all Lindblad operators is equivalent to 
showing the uniqueness of the ground state of the following positive semi-definite 
Hamiltonian (dimensionless): 

Hp = (31) 



where J" {a = ±, z or x,y,z) are the d-wave jump operators in Eq. (23). To see the 
equivalence, note that for any operator A, the matrix element (t/'IA^AI?/') > is non- 
negative. Thus, the spectrum of the Hermitian Hamiltonian is real and non-negative. 
Any zero energy eigenstate must be a ground state of the Hamiltonian. The energy is 
zero if and only if each term has zero energy. This is equivalent to J^l"^) = OVi,a. 
Hamiltonians with these properties often occur in the context of spin models, where 
they are constructed as "parent Hamiltonians" for given states [35]. We follow this 
nomenclature here. This Hamiltonian serves for the uniqueness considerations of this 
section as well as for the adiabatic passage discussed in Sec. |5| 

For later convenience, it is useful to collect the a = ± components into a 
dimensionless "reduced" parent Hamiltonian, with the normal-ordered form 

Hp = - ^{4^1,^ + c|_i^^)4_^q,_^(q+i,^ + Q_i,<,) + 2 ^ cl^Ci,^, (32) 

i,a i,(T 

where a =t, i, and the quartic terms describe effective attractive two-body interactions. 
Note that the "chemical potential" term proportional to the total particle number is 
unimportant for fixed particle number states. For a = z, the normal-ordered form 
(dimensionless) is 

i,(7 

~ X^(4+1,(T + 4-l,tT)cLcj,CT(Cj + l,<T + Cj-l.tr) + 2 ^ 4,(7^1,(7, (33) 

where the first term describes correlated hopping, and we have Hp = H^ + Hp. Following 
the discussion in the previous paragraph, we see that the d-wave state \^)d is a ground 
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state of the full parent Hamiltonian as well as of the reduced parent Hamiltonian. We 
will further demonstrate below that while the d-wave state is not the unique ground 
state of the reduced parent Hamiltonian, there are strong indications for it to be the 
unique ground state for the complete parent Hamiltonian based on symmetry arguments. 
These considerations will play an important role later in designing the implementation 
schemes. 

By construction of the jump operators, the d-wave states are ground states of the 
parent Hamiltonian. For any symmetry operation, i.e. a unitary transformation T, that 



leaves the parent Hamiltonian Eq. (31) invariant TH^T = Hp, a necessary but not 
sufficient condition for the d-wave state to be the unique ground state of the parent 
Hamiltonian is: 

r|^), = exp(«0T)|^)d, (34) 
where 0t is the phase imposed onto the d-wave state by the unitary operation T. 



The parent Hamiltonian Eq. (31) has two obvious global symmetries: phase 
rotation invariance associated with particle number conservation, and translational 
invariance associated with momentum conservation. In the following, we discuss these 
symmetries: 

Phase rotation invariance - The symmetry is generated by T<^ = exp iipN, where the 
number operator N = ^ cj^^Cj^o-- Since the d-wave state is an eigenstate for both Hp 
and N, with Hp\'$)a = and N\'^)d = 2A^, it is also an eigenstate of with eigenvalue 
exp 21^9 A^. Thus for a given fixed particle number no degeneracies occur according to 
the above criterion. 

Translation invariance - The symmetry is generated by the total center-of-mass 
momentum operator in 2-D: 



2,(T 



Tr = exp irP, (36) 

for which T^HpT'^ = Hp. It is easy to check that P\'^)d = 0, i.e. the ci-wave state 
is a momentum eigenstate with zero eigenvalue. Therefore, the translational symmetry 
does not lead to degeneracies. Note the close relation of the terms in the momentum 
operator to the jump operators, which may be transformed into each other by fiipping 
the spin on the site in the middle and changing from the symmetric to the antisymmetric 
superposition on sites {i ± 6^,1 ± Cy}. 

Discrete symmetries - On the bipartite lattices, and only for those, we find an 



additional discrete symmetry for the reduced parent Hamiltonian Eq. (32), but not 



the complete parent Hamiltonian Eq. (31). A bipartite lattice can be split into two 



equivalent sublattices {A, B) in such a way that each lattice site of A{B) is surrounded 
by lattice sites of B{A). Examples are equally spaced lattices with even number of sites 
and periodic boundary conditions in one dimension, or even site square lattices with 




Figure 4. Fidelity and entropy evolution of the master equations for 4 atoms on a 
1-D chain with 4 sites, (a) The fidelity is with respect to an antiferromagnetic Neel 
state. The dashed curve represents the evolution of the fidelity with respect to the 
other antiferromagnetic state of the system with a total spin flip; (b) The fidelity is 
with respect to a f-D singlet pairing state. The dashed curve shows the evolution 
without { Jf } jump operators. 



periodic boundary conditions in two dimensions. Given a certain finite lattice in any 
dimension, it is straightforward to clieck bipartiteness. 
On tlie bipartite lattices we find the symmetry: 

Td : Ci,t -Ci,t; Q,; Ci,; for i e A, 

Ci,t Ci,t; Ci^i Ci^i for i e B, 

and analogous for c|_^, so that T^HpT^^ = Hp, while T^HpT^^ ^ Hp. This 
transformation is canonical. A second quantised representation of the symmetry is 
highly non-local, similar to Shiba transformations for the Fermi-Hubbard model [23] , 
but in this context we only need its action on the Hamiltonian and the state. Applying 
the transformation to the rf-wave state we have 

= J](T,Dn'^kac), T,A— E'^-/(^'«)^M+e., (37) 

i V 

where /(z, a) = 1 for a = +, i G -B and for a = — , i G A\ f{i, a) = —1 for a = +,i E A 
and for a = —,i G B. Thus the (i-wave is not an eigenstate of Td which implies 
degeneracy under Hp but not Hp. The degeneracy emerging from this is two-fold for 
any lattice size, which we will see in later sections. 

In general these symmetries help to classify the lattice configurations (especially 
for finite lattices) under which we can expect uniqueness. Provided that there are no 
other symmetries of the full parent Hamiltonian that transform the d-wave state into 
other distinctive states, the d-wave state should be the unique dark state. Though 
we cannot rule out constructively the existence of other symmetries which may bring 
in additional degeneracies, the analysis here provides useful insights on the possible 
existence of degeneracy. 
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3.2. Uniqueness of the steady state: Numerical simulations 

To further verify the uniqueness of the dark state, and more importantly, the uniqueness 
of the steady state for the dissipative dynamics given by the jump operators, we have 
performed numerical simulations of the master equation dynamics on finite size systems. 
Due to the translational symmetry of the pairing states, we impose periodic boundary 
conditions on the finite 1-D and 2-D systems. We have also taken periodic boundary 
conditions on the jump operators to reduce finite size effects and to be consistent with 
the definition of the pairing states on a finite lattice. In sufficiently large systems, the 
jump operators drive the atoms in the bulk into the desired phase, and the mixing of 
states at the boundary becomes negligible so long as the system is in the thermodynamic 
limit. 

We have evolved the master equations for antiferromagnetic Neel states and singlet 
pairing states for finite size systems in 1-D. The results are shown in Fig. |4| which 
clearly demonstrate the uniqueness of the steady state anticipated above. Both the 
fidelity and entropy evolution indicates that in both cases the system is driven into the 
desired state regardless of the initial state. In comparison, when only jump operators 
of the reduced parent Hamiltonian are applied (dotted curve in Fig. gb)), the final 
fidelity approaches 0.5. For d-wave pairing states in 2-D, we have carried out a quantum 
trajectory simulation for small plaquettes. The evolution of the fidelity with respect to 
the d-wave pairing state indicates that the system approaches the final pure steady state 
exponentially, which implies the existence of a dissipative gap, analogous to the energy 
gap for the ground state of the BCS pairing state (see Fig. [5]^a)). While in a finite 
system, a gap of the order is expected due to the finite linear dimension L, in the 
thermodynamic limit, this gap vanishes. In the next section, we will derive a mean- 
field theory of the master equation in the thermodynamic limit, which shows that a 
dissipative gap appears naturally from the mean-field expansion of the master equation 
for pairing states, demonstrating the dissipative gap in our numerical simulations is not 
a mere finite size effect. Similar results can be obtained for the paired states of p-wave 
symmetry for spinless fermions (see Fig. jsj^b)). 

4. Mean-field expansion of the master equation and the dissipative gap 

In this section, we will develop a mean-field theory of the master equation for driven- 
dissipative pairing states in the thermodynamic limit which is valid at late times, where 
the system is close to the BCS-type pairing state. For this purpose, we switch from the 
fixed-number state representation to the fixed-phase (coherent state) representation. 
The justification for this procedure builds on two properties. First, the exactly known 
fixed-number pairing dark states discussed above exhibit phase locking among different 
fermion pairs. Such a property is reproduced in the coherent state representation 
\i^)d ^„ " I vac) = exp (aD^) |vac) (a = |a|e*^), with a fixed global phase 6 [c.f. 
Eq. (It])]. Second, for the BCS wavefunctions under consideration, these representations 
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Figure 5. Quantum trajectory evolution of tlie master equations for d-wave and p- 
wave states in 2 dimensions, (a) Evolution with d-wave jump operators on a 2x6 
ladder with 4 atoms; (b) evolution with p-wave jump operators on a 4x4 plaquette 
with 4 atoms. The insets indicate the existence of dissipative gaps in both cases, which 
render the convergence to the steady states exponentially fast. This result is robust 
in the thermodynamic limit as revealed by our mean-field theory. The fidelity (solid) 
is calculated by averaging over 1000 trajectories. These trajectories are then bunched 
into 100-trajectory groups, whose standard deviations are then calculated to show the 
sampling errors (dashed). 



are equivalent in the thermodynamic hmit, which can be verified exphcitly from a 
consideration of the relative number fluctuations in the flxed-phase state. Indeed, one 

1/2 

flnds for the variance AA^ = ( ,l>4^^ ) ~ -y^, where N is the total particle number 



operator, the average is taken with respect to the flxed-phase BCS state, and N is the 
number of degrees of freedom. The general normalised wavefunction of the BCS pairing 
state then reads 

1 *^ I I 

= Ilf /.^i I. + ^==it^^,JI-ae), (38) 



where (p^ is the pair wavefunction in momentum space, 6 is the overall phase of the 
pairing state, |a| is a real number flxing the average particle number density, and the 
product over q runs over the flrst Brillouin zone. Note that for 1-D singlet state, 
(Pq = cosg, while for 2-D d-wave state, ip^ = cos{qx) — cos{qy). The flxed-number d-wave 
state can be obtained from the coherent state representation via a number projection, 

l*)<^= / ^e''^\Di9)) = iDY\y^o), (39) 



27r 

where N is the total particle number corresponding to the density 

fdcij^ap^ (40) 

In the following, we will take a = 1 for simplicity, which corresponds to an initial state 
with given total particle density n = J i^^\2 - For initial states with different 
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total particle densities, one just needs to substitute in the following discussions 



with aipq, where |a| is determined from the number equation Eq. (40). Finally, we 
note that in the grand canonical ensemble, the overall phase 6 of the pairing state 
can take any value. Hence we will take ^ = in the following discussion and write 
\D) = \D{0)). Indeed, the phase 9 is not determined by the microscopic dynamics, 
since the jump operators are particle number conserving (charge 0). Therefore, the 
phase will be chosen spontaneously in the thermodynamics limit. This is the analog of 
the spontaneous symmetry breaking in the dissipative context. 

The d-wave state has three non-zero bilinear expectation values for each momentum 
mode relevant to the corresponding quantities in the subsequent mean-field theory, i.e. 
particle number and order parameter: 

I 1 2 

= (/^|4q,.C±q,.|D) = ^^j^, (41) 

= (D|c±q,tc^q,4|i?) = -nnr^' (^2) 

-L -r IV-'ql 

A; = (/^|cUA,ti^) = -rn|;p' (43) 

where a =t, i- All other expectation values vanish on the above state. 

A mean-field theory analogous to the BCS approximation for superconductivity 
can be set up based on the proximity of the density matrix to the d-wave state, giving 
rise to an ordering principle for devising a controlled mean-field approximation in the 
late time evolution, which is thus useful to study the final stages of the master equation 
evolution. 

Starting from the fact that the coherent representation of the pairing state is a 
product in momentum space, we will require this property also for an approximate 
ansatz for the density matrix, p = pq for the solution of the master equation at late 
times. This may be viewed as a Gutzwiller factorisation approach in momentum space, 
which has been used previously for a mean-field decoupling of bosonic master equations 
[T7] . This ansatz will enable us to derive a late time master equation quadratic in the 
fermion operators which contains information of the complex excitation spectrum, i.e. 
the damping of the lowest fermionic single particle excitations. 

To implement the approximation, we Fourier transform the jump operators to 
momentum space: 

■^fc" = E ¥''i<^%-k, « = ±,^, (44) 
q 

where (yjq is the pair wave function and reflects the similarity in construction of the 
jump operator and the corresponding pairing state. This gives rise to a Liouvillian in 
momentum space 

c[p] = J2 {^JM^ - {JVJ^^ P}) ■ (45) 
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The mean-field product ansatz for tlie density operator is 

P = n^q' /^q = *^^q/^' t^^Pq = l^q) (46) 
q 

wfiere eacli Pq spans tlie subspace for {q 'I', — q 1} ^7\. We tfius allow for a residual 
entanglement of the momentum modes {q, — q}, as in the BCS treatment, which is 
necessary to describe pairing. The second equation is the projection prescription to 
obtain the density operator for each mode which can be used to derive the equation of 
motion dtPq. 

We take the partial trace tr^q on both sides of the master equation, for which 
we trace over all degrees of freedom outside the subspace {q ti ~Q i}- Then, in 
the spirit of BCS theory, we choose the relevant mean fields, i.e. macroscopically 
occupied expectation values in the dark state. As mean fields of linear and trilinear 
correlations vanish due to the Fermi statistics, and that of the quartic and higher order 
correlations connect momentum modes different from ±q, and thus are small compared 
to the macroscopic expectation values close to the steady state in the thermodynamic 
limit, we keep only mean fields of quadratic correlations, i.e. density mean fields or 
condensate with zero center of mass momentum [cf. Eqs. (41) (42) (43)]. We thus 
see how the proximity to the steady state can be used as an ordering principle for a 
mean-field theory at late times, which is based on the exact knowledge of the fixed- 
phase steady state density matrix. Note that we use the commutation (as opposed to 
anti-commutation) properties during the process 

[Pp> Pq] = [Pp> Cq] = [Pp, cl] =0 for p ^ q, (47) 

which is equivalent to the assumption that there is an even number of fermions in this 
mode. 

We then have the following results: 



tr^qE(24Wk^^-{4^4^P}^ 

k ^ 
2A{1 + Iv^qH X |7q,;Pq7l4 " ^ {7L7q,; , Pq} } , 

tr^q J] (24-pJk ^ - {4"^4^P}) = 

k 

2A{1 + l^qH X |7q,tPq7q,t " ^ {7i,t^q,t ' Pq) | > 

tr^q 5^ (2 J^pJ^t _ {J^tj.^^ij ^ 2A{1 + Iv^qp) 

k 

X |7q,tPq7i,t + 7q,;Pq7i4 " ^ill^lnA + Pqlj ' 

where A = J (^^jqq^^ ^ 0, and the integration runs over the first Brillouin zone. 



(4J 
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Of particular interest and importance are the {7q,cr} operators, which coincide with 
the definition of the Bogohubov quasiparticles: 

Vl + bql 

Indeed, the mean -ield BCS pairing state is the vacuum state for the Bogohubov 
quasiparticles, 

7q,t|/^) = 7-q,tl^) = 7q,4l^) = 7-q,;l^) = 0. (50) 

Furthermore, when properly normalised as in Eq. (49), these quasiparticle operators 
obey the closed fermion algebra: 

{7q,<7,7i',a'} = Sci,ci'5„,a', 

{7q,a, 7q',a'} = {7!,^, 7^,,,'} = 0, (51) 

i.e. they are related to the original fermion operators by a canonical transformation. 
We note that the operators J^'^ do not exhibit a closed algebra structure. 



From the expressions in Eq. (48), we may identify a damping spectrum in the 
dissipative part of the master equation: 

K^ = 2AK{l + \ip^\-'), (52) 

with A = [1 — I/V2) in 1-D, and A ~ 0.36 in 2-D. It is important to note that the 
damping spectrum is gapped, i.e. /tq > 2Ak, is bounded from below. This behavior 
exhibits strong parallels to the equilibrium problem of paired fermions, where pairing 
is protected by an energy gap. Furthermore, it has the important implication that the 
approach to the d-wave dark state will be exponentially fast, in contrast to a bosonic 
system where wavelengths of arbitrary length make the approach to the dark states with 
long range order polynomially fast only [12]. This result is refiected in the quantum 
trajectory simulations in the previous section. 

Similarly, for the parent Hamiltonian of the pairing state, it is straightforward to 
derive 

tr^q[iJp, p] = /«q7i,a7q,<x, Pq] • (53) 

Clearly, the dissipative mean-field theory developed here can be applied to pairing 
states with other spatial symmetries. As an important example, we discuss the result 
for complex p-wave pairing states. The p-wave pairing state for spinless fermions can 
be written as: 

\P) = n( . ^ = + —F^^^=ciclj\va.c), (54) 

where q runs over half of first Brillouin zone, e.g. qx > 0. The pairing wavefunction 
= (/)_q = — fq, whlch Is rcqulrcd for p-wave symmetry: in 1-D, (fg = 2i sin q; in 2-D, 
(fq = 2i{smqx ± ising^). 
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Following the previous derivations, we define the momentum space jump operators: 

= ^ V5qC]jCq_k, (55) 

q 

where the summation over q is over the first Brillouin zone. Taking the partial trace 
over the degrees of freedom outside the subspace {q, — q} (q spans half the momentum 
space here, e.g. > 0), and identify the mean fields as before, we may arrive at the 
effective Hamiltonian 

i 



Hpff = — 

q 



(^i^'i + 7lq7-q) , (56) 



where the summation of q is over half the first Brillouin zone {q^ > 0). The Bogoliubov 
quasi-particles are given as 

^'l = /I 12 ^^'l ~ 2^qcLq) . (57) 

The dissipative coefficient is given as 

K^ = AK{l + \2ip^\^), (58) 

where A = J i+py p ' with the integral running over the first Brillouin zone. We 
find A ~ 0.19 in 1-D {d= 1), and A ~ 0.23 in 2-D (rf = 2). 

Finally, we note that the fermionic quasi-particle operators in Eq. (57) formally 
correspond to the Bogoliubov operators of a p-wave Hamiltonian in the limit /i — >^ — oo 
(/X is the chemical potential), which means that the state is in the strong pairing limit 
describing a state of delocalised tightly bound molecular pairs and is topologically trivial 

The generation of stable topological order is however possible in a modified setting, 
as has been established recently [20] . 



5. Adiabatic passage to the ground state of the Hubbard model 

As argued above, the pure mean-field state with the correct symmetry (antiferromag- 
netic at half-filling and d-wave otherwise) is a convenient initial state for the quantum 
simulation of the ground state of the Fermi-Hubbard model (FHM). With a suitable 
adiabatic passage, it would be possible to connect the pure dark state with the ground 
state of the FHM. A first guess as to how to reach this ground state based on the expe- 
rience with purely Hamiltonian systems might be a simple adiabatic passage, in which 
the Liouvillian is switched off while the FH Hamiltonian is ramped up. Small scale nu- 
merical simulations of the time evolution suggest that this procedure does not work for 
our combined system with both unitary and dissipative evolution. This is understood 
from the fact that the mean-field state is not an exact eigenstate of the FH Hamiltonian, 
such that unitary and dissipative evolutions compete. As a result of this competition, 
the steady state density matrix in general describes a mixed state, instead of a pure, 
zero entropy state [IE]- We thus observe that a dissipative gap cannot play the 
role of an energy gap in standard adiabatic passage schemes. We therefore propose a 
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Figure 6. Adiabatic passage connecting the antiferromagnetic state and the mean- 
field d-wave state with the ground state of the Fermi-Hubbard model, (a) The initial 
state is an antiferromagnetic Neel state on a 2 x 2 plaquette with 4 atoms; (b) the initial 
state is a d-wave state on a 2 x 4 ladder with 4 atoms. We calculate the evolution of the 
fidelity of the instantaneous system state with respect to the final ground state of the 
FHM. (inset): Time dependence of the ramping parameters h{t)/Um and J{t)/Um- 
The interaction energy U{t) in the FHM linearly increases from to its maximum 
value Urn during the ramp (not shown), with the final state corresponding to a strongly 
correlated situation with J /U = 0.1. 



modified adiabatic passage, which uses the parent Hamiltonian of the mean-field state. 



constructed from the complete set of jump operators as in Eq. (31). By construction, 
this parent Hamiltonian has the mean-field state as a gapped ground state, and therefore 
provides for an energetic stabilization. The passage now proceeds by first turning off 
the dissipation while the parent Hamiltonian is applied, then simultaneously ramping 
down the parameters of the parent Hamiltonian while ramping up the parameters of 
the FHM. In this way, as long as the symmetry patterns of the mean-field target state 
and the ground state of the FHM Hamiltonian are the same and no phase transition is 
crossed, an energy gap persists through the whole passage. Indeed, we show numerically 
that this modified adiabatic passage ensures an efficient transfer into the desired ground 
state. 

For the antiferromagnetic Neel state at half-filling, the parent Hamiltonian 
(dimensionless) reads 

^' = E(^fi'-^£+4.^M)+jIji, (59) 



where the jump operators jf^ and ji are defined in Eqs. (19, 21). We have performed 
numerical simulations of such an adiabatic passage for a 2 x 2 plaquette with 4 atoms. 
The result is shown in Fig. [6]^a). Indeed the initial Neel state can be adiabatically 
connected with the ground state of FHM at half-filling with high fidelity. 



For the case of d-wave state, the parent Hamiltonian in Eq. (31) by construction 



has the initial d-wave state as an exact eigenstate and thus supports the d-wave state 



Driven-dissipative many-body pairing states 



22 



obtained from the dissipative evolution. From Eq. (53), it is clear that the single 
fermion excitations on the d-wave state are gapped if the system is sufficiently far away 
from half-filling. As a consequence, all requirements for an efficient adiabatic passage 
are met. Note that single fermion excitations above the ground state manifold of the 
reduced parent Hamiltonian are also suppressed by an energy gap, as if^^ = ^-ffefr on 
the mean-field level. We will make use of this important fact in Sec. 6^ to design a 
modified adiabatic passage. 

The time-dependent Hamiltonian describing the adiabatic passages is given as 

Hit) = h{t)H, + U{t) cltC^t44C4 - Jit) E ^l^^- (60) 

where the time dependent coefficients h{t),U{t), J{t) give the precise path of the 
adiabatic passage. In practice, the time dependence of these coefficients is given by 
the rate at which the lattice potential giving rise to the FHM is ramped up, as well 
as by the rate of the effective interaction given by the parent Hamiltonian Hp. Here, 
for simplicity, we have chosen linear ramps for these coefficients (see insets of Fig. [6]), 
which already give a clear physical picture of the adiabatic passage. In practice these 
ramps could be further optimised, so that higher fidelities can be achieved in shorter 
ramp times. Consistent with the previous discussion, the role of the parent Hamiltonian 
is to provide an energy gap, and hence energetically stabilise the adiabatic passage. 

We have performed numerical simulations of the adiabatic passage with various 
finite size systems. To avoid degeneracies of the ground state of FHM due to finite size 
effects, we have taken open boundary conditions for the Hubbard Hamiltonian during 
the adiabatic process, while we retain the periodic boundary conditions for the definition 
of the initial pairing state and for the jump operators. We expect that the mean-field 
d-wave state should be efficiently connected to the ground state of the Fermi-Hubbard 
Hamiltonian so long as the d-wave symmetry of the ground state of the Fermi-Hubbard 
model is present and not completely destroyed by the finite size effect. We find that 
this is the case for ladder systems. A typical result of our simulation on a finite ladder 
is shown in Fig. |6](b). For systems in the thermodynamic limit, the symmetry property 
of the ground state is not affected by the boundary effect, and we expect an efficient 
adiabatic passage so long as the symmetries of the ground state are the same as the 
dissipatively driven initial state [381 139| HO]. 



6. Physical implementation and modified adiabatic passage 

As an illustrative example, we now discuss a proof-of-principle implementation of the 
single-particle jump operators. The scheme we describe in this section is stroboscopic, 
and involves realising the action of the jump operators in a series of steps. Though 
non-trivial to implement in present experiments, this example is made up of elements 
that are presently accessible in experiments. The example illustrates how the properties 
of the operators appearing in the previous sections, specifically that they are quasi- 
local, conserve particle number, and can be implemented based on single-particle 
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Figure 7. Level scheme using ^^^Yb atoms. The physical spin state is encoded in 
the nuclear spin sublevels of the ^Sq manifold. The spin flip operation is implemented 
via off-resonant coherent coupling to the '^Pq manifold with circularly polarised light 
(red arrows). The long lived '^Pq states are coupled to the ^Pi level in a two-photon 
process, from which spontaneous emission into a cavity is induced, leading back to the 
^So manifold. 

operations, make them favourable for experimental implementation. For an alternative 
non-stroboscopic, i.e. "always-on" continuous implementation, which is applicable in 
the case of spinless (spin-polarised) fermions such as the p-wave case discussed above, 
see 120]. 

Our example takes advantage of the properties of alkaline-earth-like atoms [HI |121 
SSI HU SSI US] . With two valence electrons, these atoms possess metastable triplet levels, 
and fermionic isotopes have non-zero nuclear spin (e.g., J = 1/2 for ^^^Yb, which we 
will choose here). This nuclear spin acts as an independent degree of freedom in the 
ground ^Sq and lowest excited ^Pq manifolds. Here, the nuclear spin will play the role 
of the physical fermionic spin degree of freedom, and the ^Pq manifold will be used as 
an intermediate state in the dissipative process. These states are depicted in Fig. [7) 
Note that as ^Sq and ^Pq are optically separated, they can be trapped in independent 
lattices using dipole traps at different wavelengths |17]. 

As a simple example, we will first discuss the implementation of jump operators 
for driving the system into the antiferromagnetic Neel state at half-filling. We will then 
move on to the more complicated cases of pairing states. 

6.1. Antiferromagnetic Neel state 

Since the Lindblad jump operators for the Neel states act on unit cells of two sites, 
we will focus on operations on two adjacent sites i and j. These will be carried out in 
parallel on pairs of adjacent sites along the lattice. As illustrated in Fig. [7[ spins states 
are encoded in the nuclear spin sublevels of the ^Sq manifold. We further assume that 
atoms in the ^Sq and ^Pq manifolds are trapped in independent optical lattices, and 
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(a) Step(l) Step (2-3) Step (4) 




Figure 8. Implementation sctieme for tlie jump operator cj ^c^,^ for 1-D 
antiferromagnetic Neel state, (a) Tiie decay cliannel given by the jump operator is 
not blocked; (b) the decay channel is Pauli blocked. 

that ^Po is trapped in a superlattice with a period of two sites, defining pairs of sites, 
where we label the left well i and the right well j. Initially, the superlattice potential 
is arranged in such a way that the potential well at site i is much deeper than at site j 
(with an energy difference of the lowest state in each well of the order of several kHz). 
The action of the jump operator c'j-^Ci^i can then be realised by performing the following 
operations (see Fig. |8]): (1) apply a circularly polarised vr-pulse selectively on site i 
coupling the ^Sq and '^Pq manifolds so that any atom originally in the state | l,^ So) 
will end up in | '|','^Po); (2) adiabatically manipulate the superlattice potential so that 
the population at site i is transferred to site j; (3) couple the ^Pq and ^Pi manifolds 
off-resonantly, so that the population in |'|',^ Pq) should decay to 1 1/ So) on site j if and 
only if the state on site j is empty; (4) repeat the steps (2) and then (1) (in reverse 
order) to bring any remaining population in ^Pq back to the ground state manifold. 

Before moving on to extend the scheme to jump operators associated with 
pairing states, several comments are in order: (a) the jump operator is implemented 
stroboscopically, which places requirements on the time scale of each step of operations 
listed above, such that the total time of evolution should be much longer than the time 
scale of operations; (b) the jump operators are implemented in parallel for each pair of 
lattice sites along the lattice; (c) during the excitation of the population from ^Sq to 
^Po, as the line-width of the metastable '^Po is on the order of lOmHz for ^''^Yb, the 
bias between different subwells in the superlattice potential ensures site selectivity; (d) 
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Figure 9. Implementation scheme for jump operators of 1-D d-wave state, (a) Only 
one of the decay channels is blocked; (b) both decay channels are Pauli blocked, 
therefore the jump operator does not change the system configuration. The state 
is a local dark state for this local jump operator. 

the nuclear spin is conserved during the decay process, which is guaranteed by the large 
detuning from the ^Pi manifold. The nuclear spin conservation can also be realised in 
this case by applying a large magnetic field so that electronic spin and nuclear spin are 
decoupled [Hj. 

6.2. d-wave pairing state 

We now extend the ideas above to the implementation of jump operators for driving 
the system into d-wave pairing states. For the d-wave jump operators, an additional 
constraint on the dissipative process is that atoms on quasi-local sites, e.g. site i -\- Cx 
and i — Cx, should decay coherently. To satisfy this requirement, we couple the system 
to a cavity with a finite linewidth. An atom (or atoms) at the sites i + e^ and i — e^ will 
then be coupled collectively to the cavity mode, ensuring that the decay is coherent. For 
clarity, we first describe the implementation procedures in 1-D, and choose the example 
of = (c|_,_i^ + c|_^^)ci,4.. The step-by-step implementation scheme is shown in Fig. 
[oj (1) We first assume that the '^Pq state is initially trapped in a lattice of three times 
the period as that for the ^Sq state, defining blocks of three sites in the original lattice. 
Using this, we excite any spin-down atom in ^Sq on central site to the spin- up state of 
the ^Pq manifold, using cr"*" light. (2) We then add an additional potential, splitting 
this site into two, and separate these sites so that the mode of atoms confined in them 
overlap the right and left sites of the original three-site block. (3) We induce decay by 
coupling atoms in the ^Pq state off-resonantly to the ^Pi state, with coupling strength 
Q, and detuning A. If we couple the ^Sq-^Pi transition to a cavity mode with linewidth 
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Figure 10. Generalisation of the implementation of d-wave jump operators to a 2- 
D lattice. Only the manipulation of the upper superlattice is shown here, (upper 
panel) After the population in the central site of the 3-by-3 plaquette in the ^Sq level 
is excited to the superlattice of ^Pq, the potential at the central site is adiabatically 
lowered so that the state is adiabatically connected to the one where the relative phase 
between the central site and its neighbors is negative; (lower panel) the superlattice 
is then shifted adiabatically in the y direction, splitting the remaining population in 
the central site along y so that the correct relative phase with d-wave symmetry is 



imposed as given by in Eq. ( 22 ) . One may then follow the procedure for 1-D singlet 
pairing state implementation. 



r and vacuum Rabi frequency g, then the decay will be coherent over the triple of sites. 
In the limit A ^ r2,r and T ^ max(^, ^, ^), we adiabatically eliminate the cavity 
mode and the intermediate far off-resonant sta te ^Pi, and ob tain an effective decay rate 
Feff = ~ 9kHz for typical parameters (see Appendix C). Note that Fermi statistics 
will be observed in this process, and that we assume that the atoms remain in the 
lowest band, as all parameters are smaller than the trapping frequency in the lattice 
(see Ref. [31] for more details of Pauli-blocking of spontaneous emissions in this sense). 

Other jump operators, J~ and J^ can be implemented by applying rotations in the 
nuclear spin before and after the three steps above. For J^~, one exchanges the spins 
with a TT-pulse, whereas for Jf , one must apply a 7r/4 rotation in the nuclear spin basis 
before and after the operation. In addition, for J/, both spin states should be excited, 
and coherence of nuclear spins is maintained throughout the operation. This can be 
achieved by either going far off-resonant for the field coupling '^Po and ^Pi manifold, or 
by applying a large magnetic field as described in Ref. [H]. 

This scheme can be generalised to 2-D by considering 3-by-3 plaquettes defined by 
the appropriate superlattice potential for the ^Pq level. As in the 1-D case, we require 
an adiabatic manipulation of this potential in step (ii), although here the depths of the 
wells must be adjusted to ensure that the correct relative phases are obtained for atoms 
"transported" in different directions (see Fig. 10 and its caption). 
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Figure 11. Implementation of {jf )'^ jf in tlie parent Hamiltonian for 1-D case. 
Firstly, the population of the outer sites are excited to the upper lattice. The super 
lattice is then adiabatically tuned from a double well structure to a single well, so that 
the state \L) + \R) is projected to the lowest level of the single well potential. The 
interaction is then induced via a Feshbach resonance for instance after exciting the 
population of the opposite spin in the central site to the superlattice. 



6. 3. Implementing the reduced parent Hamiltonian and modified adiabatic passage 

Here, we extend the scheme above to implement the reduced parent Hamiltonian for the 
d-wave state stroboscopically. We see from the discussion in Sec. |3]that the mean-field 
d-wave state is in the ground state manifold of the reduced parent Hamiltonian. As 
we have discussed in Sec. [sj this degenerate ground state manifold (two-fold in 1-D, 
four-fold in 2-D) is protected by an energy gap from single fermion excitations under the 
reduced parent Hamiltonian. Furthermore, we will also show below that an adiabatic 
passage with high fidelity can be achieved by a modified adiabatic passage scheme. 

The implementation of the reduced parent Hamiltonian is similar to that of the 
jump operators, except that the dissipative part is replaced by an induced phase shift. 



As shown in Eq. (32), the reduced parent Hamiltonian contains an effective interaction 
term and a term proportional to the total particle number that is not important for states 
with fixed particle number. To implement the effective interaction term stroboscopically, 



as illustrated in Fig. 11 with the example of 1-D Hamiltonian {J^yj^, the following 
steps are required: (1) any spin down atoms in the left and right well of the ground 
state lattice potential are transferred to the superlattice potential of the ^Pq manifold; 
(2) the double-well in the superlattice potential is merged into a single well, during 
which process the symmetric state in the double-well potential is mapped to the lowest 
motional state of the final single well potential; (3) a phase shift is then induced to 
generate an interaction only if the spin-down state in the ^Sq manifold and the spin-up 
state in the ^Pq manifold are simultaneously occupied. This can be achieved, e.g., by 
applying a tt pulse between the spin-down state in the ^Sq state and the spin-down state 
in the lowest motional state in the superlattice potential, and inducing an interaction 
between the different spin states in the ^Pq manifold via an optical Feshbach resonance. 
Finally, to implement {J^YJ^, the spins should be exchanged while the above procedure 
is carried out. 
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Figure 12. Adiabatic passage connecting the mean-field d-wave state with the ground 
state of the Fermi-Hubbard model on a 2 x 6 ladder with 4 atoms with (a) complete 
parent Hamiltonian; (b) reduced parent Hamiltonian. We calculate the evolution of 
fidelity of the system state with respect to the final ground state of the FHM. (inset): 
Time dependence of the interaction rate U{t)/Um- In both cases, h{t)/Um is ramped 
down linearly from 0.05 to 0, while J{t)/Um is ramped up from to 0.1. Similar to 
Fig. [6]^b), the fidelity remains small until late in the adiabatic process. The overlap 
with the ground state of the FHM only becomes large at late times when the repulsive 
interaction is large enough to overcome double-occupancies that occur due to the form 
of the mean-field pairing wavef unction. 



With only the reduced parent Hamiltonian, we find that given an optimised ramping 
scheme, the mean-field d-wave state can still be adiabatically connected with the ground 
state of FHM on small lattices. Fig. 12 shows such an example, where ramps with the 
complete parent Hamiltonian and with the reduced parent Hamiltonian are numerically 
simulated for 4 atoms on a 2x6 ladder. For the adiabatic passage with the reduced 
parent Hamiltonian, we ramp up J{t) and U{t) separately. In both cases, we have very 
high fidelity at the end of the ramp. For ramps with the reduced parent Hamiltonian, 
the high fidelity is due to the large overlap (~ 0.95 for most ladder systems) between 
the mean-field d-wave state and the ground state of the time dependent Hamiltonian 
at the beginning of the ramping process (when there is no on-site interaction). For the 
numerical simulations that we considered here, this overlap also sets the upper bound 
for the final fidelity of the ramps. 



7. Conclusions 



We have proposed an approach for the preparation of many-body pairing states of given 
symmetry for fermionic atoms in an optical lattice via driven dissipative processes based 
on suitable reservoir engineering. We have discussed in detail the strategy of designing 
the jump operators making use of the Fermi statistics, which gives rise to the dissipative 
preparation of the initial pairing state. This process is in general efficient, due to the 
existence of a dissipative gap. We then argued for the uniqueness of the pairing state as 
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the steady state of the dissipative dynamics, both from symmetry considerations, and 
via small scale numerical simulations. Note that for realistic finite size systems such 
as plaquette geometries |18], it is also possible to design jump operators for specific 
many-body states defined on the finite system, in which case one may need to design 
special "boundary" jump operators to make the state unique. 

We then discussed the adiabatic passage process that could be used to connect 
the driven-dissipative mean-field state with the ground state of the FHM. As our d- 
wave state is not an eigenstate of the FHM, directly ramping down the dissipation 
rate while ramping up the FHM leads to competition between the coherent and 
dissipative dynamics which would not drive the system into the ground state. We 
therefore introduced the parent Hamiltonian of the d-wave state, a semi-positive 
Hermitian Hamiltonian constructed from the jump operators. By construction, the 
parent Hamiltonian has the dark state of the dissipative process as its ground state. We 
illustrated via small scale numerical simulations that the ground state of the FHM can be 
adiabatically connected with the mean-field state of the relevant symmetry via optimised 
adiabatic paths. This is in similar spirit to the recent experimental demonstration of 
antiferromagnetic order in an optical lattice [S] , where the desired eigenstate of the Ising 
model is prepared via adiabatic passage from a starting state that has low entropy and 
sufficient overlap with the final state. We note that it is possible to extend these small 
scale numerical calculations by applying time dependent density matrix renormalisation 
group (t-DMRG) methods. In fact, quantum trajectories methods could be combined 
with t-DMRG methods [H] in order to perform larger-scale simulations of the dissipative 
preparation process and the adiabatic ramp together. 

Finally, we discussed a proof-of-principle physical implementation of both the 
jump operators and the parent Hamiltonian using alkaline-earth-like atoms, which 
illustrated that the properties of the jump operators discussed here are favourable for 
implementation. We mainly focused on the implementation of d-wave pairing state, but 
similar implementations can be readily found for pairing states of other symmetries, so 
long as the jump operators are quasi-local and involve operations manipulating only 
single particles. 
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Appendix 

Appendix A. Two-particle jump operators for d-wave pairing state 

In this appendix, we derive in detail the two-particle jump operators for pairing states 
with d-wave symmetry as appeared in Eqs. ( I6p!7 ). 



Appendix A.l. One- dimensional case 

As an example, we choose 

ii = Ci+uQt - CilCi+n- (A.l) 
The commutation relations then give 



+ ('^Lit'^4 4t4-i4.) + ('^Lit'^i+24 4+2t'^I-u)| • (A.3) 
It is straightforward to show that if we define 



Xi 



('^l+it'^I+24, + 4+2t4+u) + ('^l+it'^4 ~^ ^l-t^l+u) 



+ (4-it'^4 + (4-it4+24. + 4+2t'^l-i4.) J ' (A.4) 

then XiBi = 0. 

The symmetry in the expressions above suggest that we may simplify these 
operators by choosing = Cj+^Cjf . The commutation relations then have the form: 

Ai = 1 — cl^Ci-^ — c|_^^Q+i4 — c|_^2t'^*t ~ cj-u^i+u (A. 5) 

B,= - 2(4^ + 42t)(4+u + 4-u)- (A.6) 
The most straightforward choice of Xi would be Xi = B^, as Bf = 0. More generally, 
XiBi = is satisfied so long as the pair operators in Xi can be factored out to contain 
either (cj^^ + c|_^2t) ('^I+u ~'~ This actually allows some freedom in choosing the 

remaining part of the Xi operator. 



However, in this second scenario, the existence of a constant term in Eq. (A. 5) 



renders Eq. (14) not equal to zero even if Eq. (15) is satisfied. The resulting jump 
operator would then not give the desired dark state. To solve this problem, one needs 
to introduce appropriate symmetry into the design of the jump operator. Notice that 
assuming the translational symmetry, the creation operator of the state can also be 
written as: 

Correspondingly, we examine the following factorised 

6 = (q+u - Ci-u)Qf (A.8) 
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Note that the choice of the negative sign here is to ensure that no constant terms appear 
in the expression for Ai. 

For the commutation relations, we now have 



(A.9) 
(A. 10) 



This implies that we may satisfy the dark state requirement by choosing a jump operator 
of the form 



J^ = C^MCi^, M = (4+1 , + , ) - C^^ll) 



(A.ll) 



as given in Eq. (16). 



Appendix A. 2. Two-dimensional case 
We define 

6 = (Cj+e^4, + Ci+eyl)Ci-\, 

whose commutation relations are: 
Ai = 



(A.12) 



+ (cl+2eyt'^«t " 4+2e^t'^n) + i4+e^-eyt'^i1- " 4+ey~e^t^it^ 
Bi = 2(c|+g^_g^^ — c]'_,_e^_e^^ — c|+2eyt 4+2e^t) 



Following the same derivation as in the previous section, we find 



(A.13) 
(A. 14) 
(A.15) 



where is an arbitrary superposition of single-fermion creation operators. Note that 
this is the most straightforward choice to satisfy Xi^i = 0, other solutions may still 
exist. 

Finally, we see that in the case of a d-wave state on a 2-D lattice, the jump operator 
takes the form: 



M = - pA+e^^Ci+e^i + Ci+eyi), 



(A.16) 



where p±x = 1, p±y = — 1, as given in Eq. (17) 
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Appendix B. Jump operators for the fixed-phase state 

hi this Appendix, we discuss the general formahsm for the construction of jump 
operators for the fixed-phase state, starting from number-conserving jump operators 
with known unique dark state with fixe particle number. These jump operators describe 
dissipative processes for which the total particle number is not exactly conserved, 
whereas the average particle number approaches the steady state value determined by 
the parameters of the dissipative process. In the following, we will first discuss the 
pairing states of spinful fermions, before extending the formalism to spinless fermions. 



Appendix B.l. Pairing states with spins 



We only consider separable pairing states, i.e. pairing states whose spin degrees of 
freedom can be factorised. Then the general number-conserving pairing state for spinful 
fermions can be written as: 

N 



^C/AM |vac> 

. I \ u / \ /J. 



N 



I vac). 



(B.l) 



where Cj = J2u Pt^^l+e^,^ M = -^M^+e^.^ translation invariant. Fourier 
transform the pairing operators into the momentum space. 



,4' 



(B.2) 
(B.3) 

With these, the fixed-phase 



where /k = Pue and g]^ = A^e 

correspondence of Eq. (B.l) can be written in the form of a coherent state: 

(e. c!4 



^E 



w. 



I vac) 



A/'exp ^^CjAl^ I vac) 

A/" JJ (^1 + a/kfl'-kCk^-^clk,^) I vac). 



(B.4) 



where a is a complex number carrying the phase of the pairing state. A/" is the 
normalization factor, and q runs over the first Brillouin zone in the product. Without 
loss of generality, this coherent state can be re-arranged into the standard BCS-type 
mean-field wave function: 



(B.5) 
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where the coefficients 



Vl+l"/k9-kP"' y/l+Kkg-kP" 

Apparently, the coherent state Eq. (B.5) is the vacuum for the Bogohubov 
quasiparticle operators: 



7k,t = WkCk.t - VkC-k,i, 
J 



7k,; = MkC-k,; + WkCk^^, 
as it is easy to verify the following relations, 7k,i- 



(B.6) 
(B.7) 



; = iKimc = 7-k,tl^)c = 

0. Therefore, these Bogoliubov quasiparticle operators are the momentum 



space jump operators for the fixed-phase state Eq. (B.5). Based on Eqs. (B.6, B.7), it 
is easy to find a more general form of the momentum space jump operators 



7k,t 



^^k t 
— r 



Mk 



'k,t 



(B.8) 
(B.9) 



where 93^(k) are arbitrary functions of k. Fourier transforming Eqs. (B.8 B.9) back to 
the coordinate space, we immediately get the quasi-local jump operators that we look 
for. 

As an illustrating example, let us investigate the simple case with (/k = 1, 
(^^(k) = 1, which implies the structure of the pairing state should be completely encoded 
in the spin-up degrees of freedom. The coherent state in this case becomes 



A/" JJ (^1 + a/ki-f-c^k,;) I vac). 



(B.IO) 



The corresponding momentum space jump operators are: 
7k,t = Ck,t - a/kC^k,4' 

7-k,4 = Ck,4 + «/-kcLk,f (B.12) 
We then Fourier transform the momentum space jump operators to the coordinate space. 



For the d-wave pairing state, p±^ 



-P±y 



(B.14) 

1, and we recover the fixed-phase 



jump operators for the d-wave pairing state in Sec. 2^ As a consistency check, we 
demonstrate below that the coherent state in Eq. (B.4) is a dark state of the jump 
operator in Eq. ( |B.14[ ): 



a 



a' 



I vac) 



I vac) - J] 



a 
n\ 



I vac) 
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a 



Pucl 



a 



vac) 



= 0. (B.15) 

Finally, we note that the fixed-phase jump operators presented here are 
superpositions of a creation operator and an annihilation operator, and hence do not 
conserve the total particle number. The average particle number on the other hand, 
is driven to the final steady state value during the dissipative process. The average 
particle number in the final steady state is given as: 



k,(T k 



(B.16) 



Appendix B.2. Spinless pairing states 



One can easily extend the formalism above to the pairing states of spinless fermions, 
with modifications to the pairing parameters due to the triplet pairing symmetry. The 
pairing state with fixed phase can be written as 



A/'exp 



a 



EE 



PuC. 




I vac) 



J]^' (^Uk + t^kCkC^k) |vac). 



(B.17) 



where indicates that k only runs over half of the Brillouin zone, e.g. > 0. 
The coefficients here are similar to the spinful case, with Uk 



^^k 



2"/k9-k 



^l+|2a/kg_kP 



/k = E,.P<^e ''^■^•^ and 



^k 



v/l+|2Q/k<7-k|2 

^ A^e^*'^ '''' . As a simple example, we 



consider the case with (7„k = 1- The coherent state then becomes 
\i^)p = ^fll'[l + 2a^clcl^) |vac). 



(B.18) 



where due to the triplet pairing symmetry, /_k = /k = ~/k- 

Following the previous approach, the general form of the momentum space jump 
operators are 

7k = ¥^(k)(ck-2«/kclj, (B.19) 

where v'(k) is an arbitrary function of k. For simplicity, we take '/'(k) = 1, and the 
Fourier transform of Eq. ( |B.19[ ) gives the quasi- local jump operators: 

^i = c, + 2aJ2p4+..- (B.20) 



Similar to Eq. (B.15), it is straightforward to check that 7^ 
Finally, the average particle number is given by 

|2«/kP 



0. 



|2a/k 



k k 

where the summation runs over the entire first Brillouin zone. With 
—ipy = ip-y = 1, we recover the results for p-wave pairing states in Sec. 2.3 



(B.21) 

-P-x = 
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Appendix C. Derivation of the effective decay rate 

In this Appendix, we derive the effective decay rate from ^Pq to ^Sq manifold via 
cavity mode during the implementation scheme illustrated in Sec. VI. The level scheme 
and various parameters are shown in Fig. ([T]). We basically need to adiabatically 
eliminate the intermediate states in the ^Pi manifolds as well as the cavity photon 
mode sequentially. 

For the clarity of discussion, we denote \a) = |^Po,0), \b) = |^Pi,0),|c) = |^So,l), 
\d) = |1q,0), where the second index indicates the number of photons in the cavity 
mode, and we have neglected the indices for the hyperfine spins for simplicity. The 
master equation for the density matrix is then 

p = —i[H,p] — —{a'^ap + pa^a) + Tapa\ (C.l) 

where a(a^) is annihilation (creation) operator for the cavity photon mode. The 
Hamiltonian under the rotating wave approximation and appropriate rotating frame 
reads: 

H = Aalab + (^alofe + h.c.) + {gala^ + h.c), (C.2) 

where aj(i=a,b,c,d) is the annihilation operator for the corresponding level, VL is the 
effective Rabi-frequency between \a) and and g is the coupling rate between the 
cavity and the atom. 

The equations of motion become: 

(C.3) 

(C.4) 
(C.5) 
(C.6) 

(C.7) 
(C.8) 

\pbc (C.9) 

(C.IO) 
(C.ll) 
(C.12) 

where pij = {i\p\j). Physically, the detuning from ^Pi manifold should be much larger 
than the its bandwidth to avoid large spontaneous emission, therefore A ^ 28MHz 
becomes the largest energy scale in the equations. We may then adiabatically eliminate 



Paa 


= -l-{pba-pab) 






Pbb 


= - l^ipab- Pba) 


- igipcb 


-pbc) 


Pec 


= - ig{pbc - Pcb) - 






Pdd 




n 

2 \Paa 




Pab 


= iApab + igpac + i 


pbb) 


Pac 


Q 

= - ^^Pbe + igPab 


F 

~ ~^Pac 




Pbc 


= - tApbc - I— Pac 


- igipcc 


-pbb) 


Ped 


F 

= - WPbd - -^Pcd 






Pbd 


= - iApbd - igpcd ■ 


Q 

- i^Pad 




pad 


Q 

= -i^Pbd 
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^Pi manifold first, wfiidi amounts to setting — Pbc — Pba — 0. Tfie resulting equations 
of motion become: 

Paa= -i^{Pac-Pca) (C.13) 

Pec = ^^(Pac - Pea) " ^Pce (C.14) 
Pdd = Tpcc (C.15) 

Pac = ^^(Pcc - Paa) + " ^)Pac " 2P«c (C.16) 

Pad = «^Pc<i + (C.17) 

c/2 r , 

Pcd = i-^Ped + ^^Pad - 2^^' {C.18) 

where we have assumed A ^ F, and neglected terms on the order of 

We may then adiabatically eliminate the cavity photon mode by pae — Pdc — 0. 
This way, we arrive at the final equations of motion 

Pad = - ^Pad + ij^Pad (C.19) 
Paa = - TcflfPaa (C.20) 
Pdd = TcSPdd (C.21) 

where the effective decay rate T^s = and we have assumed T ^ max(^, ^, 

For typical experimental parameters: A ~ lOOMHz, Q ~ lOMHz, k ~ lOMHz, 
g ~ 3MHz, we obtain an effective decay rate Feg ~ 9kHz. 
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